Федеральная служба государственной статистики ежемесячно оценивает среднюю номинальную заработную плату в России; Известны значения с января 1993 по январь 2016. Необходимо построить прогноз на следующие три года.

Попробуем поделить на число дней в месяце: Ряд не стал более регулярным, так что вернёмся к исходным данным.

STL-декомпозиция ряда:

Оптимальное преобразование Бокса-Кокса и результат его применения:

В данном случае преобразование имеет смысл использовать, так как оно хорошо стабилизирует дисперсию. Попробуем округлить параметр и взять \(\lambda=0\): Результат практически такой же. Далее будем использовать \(\lambda=0\).

Прогноз ETS

## ETS(A,A,A) 
## 
## Call:
##  ets(y = tSeries, lambda = LambdaOpt) 
## 
##   Box-Cox transformation: lambda= 0 
## 
##   Smoothing parameters:
##     alpha = 0.8968 
##     beta  = 0.0399 
##     gamma = 0.1032 
## 
##   Initial states:
##     l = 4.6769 
##     b = 0.0031 
##     s=0.1716 0.0076 -0.0044 0.0062 -0.0011 0.008
##            0.0304 -0.0293 -0.034 -0.0173 -0.0661 -0.0716
## 
##   sigma:  0.0335
## 
##       AIC      AICc       BIC 
## -289.4724 -287.1094 -227.8641

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                        ME      RMSE       MAE        MPE     MAPE
## Training set   0.02705718  4.218196  2.568493  0.0348505 2.145357
## Test set     -20.76670383 30.891186 23.044210 -9.0295990 9.886351
##                   MASE        ACF1 Theil's U
## Training set 0.1821666 -0.06878022        NA
## Test set     1.6343766  0.89928841   1.12668

Остатки:

Достигаемые уровни значимости критерия Льюнга-Бокса для них:

Остатки коррелированы, по всей видимости, модель недостаточно хороша.

Q-Q plot и гистограмма для остатков:

Распределение имеет длинные хвосты и выброс слева.

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 2.120423110^{-16}
Несмещённость Уилкоксона не отвергается 0.7450093
Стационарность KPSS не отвергается 0.1

ARIMA

Ручной подбор модели

Исходный ряд нестационарен (p<0.01, критерий KPSS); сделаем сезонное дифференцирование: Ряд всё ещё нестационарен (p<0.01, критерий KPSS). Проведём ещё одно дифференцирование: Для полученного ряда гипотеза стационарности не отвергается (p>0.1)

Посмотрим на ACF и PACF полученного продифференцированного ряда:

На ACF значимы лаги 1 и 12, на PACF — 1, 12 и 24. Будем искать модель, оптимальную по AICc, в окрестности ARIMA(1,1,1)(2,1,1)\(_{12}\):

Модель AICc
ARIMA(1,1,1)(2,1,1)\(_{12}\) -1044.6245289
ARIMA(0,1,1)(2,1,1)\(_{12}\) -1040.8500956
ARIMA(1,1,0)(2,1,1)\(_{12}\) -1042.6126473
ARIMA(1,1,1)(2,1,0)\(_{12}\) -1027.8636046
ARIMA(1,1,1)(1,1,1)\(_{12}\) -1046.7032117
ARIMA(1,1,2)(2,1,1)\(_{12}\) -1042.5910157
ARIMA(1,1,1)(2,1,2)\(_{12}\) -1043.2935219
ARIMA(1,1,1)(3,1,1)\(_{12}\) -1045.2735494
ARIMA(0,1,0)(2,1,1)\(_{12}\) -1033.1606073
ARIMA(1,1,1)(1,1,0)\(_{12}\) -1015.0624859
ARIMA(1,1,0)(1,1,1)\(_{12}\) -1044.6899974
ARIMA(0,1,1)(1,1,1)\(_{12}\) -1042.9255562
ARIMA(1,1,1)(0,1,1)\(_{12}\) -1048.3684574
ARIMA(2,1,1)(0,1,1)\(_{12}\) -1043.7285109
ARIMA(1,1,2)(0,1,1)\(_{12}\) -1042.381999
ARIMA(1,1,1)(0,1,2)\(_{12}\) -1046.6988068

Наилучшая по AICc модель — ARIMA(1,1,1)(0,1,1)\(_{12}\). Посмотрим на её остатки: Видно, что в начале ряда остатки не определены, что логично, поскольку модель сезонная. Отрежем начало ряда остатков и проанализируем их:

Достигаемые уровни значимости критерия Льюнга-Бокса для остатков:

Q-Q plot и гистограмма:

У распределения остатков достаточно тяжёлые хвосты; скорее всего, гипотеза нормальности будет отклонена.

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 2.605680310^{-17}
Несмещённость Уилкоксона не отвергается 0.3276633
Стационарность KPSS не отвергается 0.1

Из-за отсутствия нормальности в итоговом прогнозе будем использовать предсказательные интервалы, полученные с помощью симуляции.

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                        ME      RMSE       MAE         MPE     MAPE
## Training set   0.07955178  3.716015  2.369193  0.06590879 2.013870
## Test set     -19.42347092 29.977272 21.039086 -8.35565908 9.002061
##                   MASE       ACF1 Theil's U
## Training set 0.1680315 -0.1637931        NA
## Test set     1.4921661  0.9119375  1.096755

Автоматический подбор модели

Применим функцию auto.arima:

## Series: tSeries 
## ARIMA(1,1,1)(2,1,2)[12]                    
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ar1      ma1     sar1    sar2     sma1     sma2
##       0.7348  -0.5706  -0.6775  0.0262  -0.0534  -0.4508
## s.e.  0.1377   0.1658   0.3990  0.1026   0.3952   0.2579
## 
## sigma^2 estimated as 0.001057:  log likelihood=528.87
## AIC=-1043.73   AICc=-1043.29   BIC=-1018.7

Предлагается модель ARIMA(1,1,1)(2,1,2)\(_{12}\). Её AICc выше, чем у модели, подобранной вручную. Посмотрим на её остатки:

Отрежем первые 13 отсчётов и продолжим анализ:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 2.837000110^{-17}
Несмещённость Уилкоксона не отвергается 0.389983
Стационарность KPSS не отвергается 0.1

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                        ME      RMSE       MAE         MPE     MAPE
## Training set   0.08046572  3.703396  2.358651  0.06791549 2.004611
## Test set     -20.46838493 30.868728 21.933648 -8.78402003 9.374485
##                   MASE       ACF1 Theil's U
## Training set 0.1672838 -0.1581542        NA
## Test set     1.5556116  0.9097194  1.128954

Сравним остатки двух версий аримы, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены: 

## 
##  Diebold-Mariano Test
## 
## data:  resres.auto
## DM = 0.36582, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.7148
## alternative hypothesis: two.sided

Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов.

В целом автоматическая модель сложнее, её остатки не лучше, AICc — выше, а ошибка на тесте больше, так что остановимся на модели, подобранной вручную.

Итоговое сравнение

Сравним остатки лучших моделей ARIMA и ETS, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:

## 
##  Diebold-Mariano Test
## 
## data:  resres.ets
## DM = -2.5201, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.01232
## alternative hypothesis: two.sided
## 
##  Diebold-Mariano Test
## 
## data:  resres.ets
## DM = -2.5201, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.006161
## alternative hypothesis: less

Согласно критерию Диболда-Мариано, прогнозы метода ETS значимо менее точные, поэтому в качестве финального выберем прогноз найденной вручную аримы.

Финальный прогноз

Поскольку остатки ненормальны, предсказательные интервалы строятся бутстрепом

##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Feb 2016       207.4896 202.2358 213.5791 194.05194 221.2860
## Mar 2016       219.3715 210.2357 230.3308 200.43912 240.0972
## Apr 2016       225.1880 212.1612 240.7841 199.93986 252.4175
## May 2016       223.6507 207.2247 243.2913 193.74074 257.2214
## Jun 2016       232.8283 212.8136 257.3899 195.21792 273.3823
## Jul 2016       224.2908 202.1878 251.8536 182.74209 269.4255
## Aug 2016       216.0811 192.7082 246.1042 170.48145 265.2864
## Sep 2016       218.4019 192.8091 252.1885 166.18473 271.5579
## Oct 2016       220.2080 192.1521 258.0965 161.61769 279.3173
## Nov 2016       220.9240 190.7413 261.2000 158.33086 286.4834
## Dec 2016       284.4404 241.3847 339.8941 201.62677 375.5025
## Jan 2017       208.7604 174.9442 251.5439 143.93969 279.2856
## Feb 2017       207.3798 171.9648 253.8228 140.65478 284.0678
## Mar 2017       219.3534 179.4991 271.5803 145.25588 305.4075
## Apr 2017       225.2421 181.4571 283.3184 146.59201 320.2336
## May 2017       223.7567 178.1807 285.9210 141.75326 324.4977
## Jun 2017       232.9778 183.0249 301.2620 145.22601 344.1172
## Jul 2017       224.4622 174.0226 294.1353 138.00068 338.1499
## Aug 2017       216.2652 166.0780 287.5420 130.64608 328.9515
## Sep 2017       218.6019 166.3323 294.5800 129.86302 339.3105
## Oct 2017       220.4198 164.9945 301.0726 128.74046 348.2100
## Nov 2017       221.1438 163.6116 305.1979 126.97687 355.7924
## Dec 2017       284.7302 207.1066 398.3282 162.21318 465.0166
## Jan 2018       208.9767 150.1895 295.4684 117.20492 345.3028
## Feb 2018       207.5973 148.1065 298.2987 114.57097 350.5245
## Mar 2018       219.5854 155.0142 320.9565 118.87609 378.3988
## Apr 2018       225.4818 157.3582 332.5987 118.83430 397.0380
## May 2018       223.9959 155.1545 336.4681 115.49945 403.6425
## Jun 2018       233.2277 159.0596 354.1100 118.20198 427.0991
## Jul 2018       224.7035 151.5078 346.1786 112.09896 419.3710
## Aug 2018       216.4981 144.2518 337.7193 106.89095 410.0909
## Sep 2018       218.8376 145.3020 346.4828 106.51824 420.5966
## Oct 2018       220.6577 144.6078 352.3042 105.47649 433.9250
## Nov 2018       221.3825 144.2202 357.8302 103.73521 449.3872
## Dec 2018       285.0378 183.8538 466.8318 130.25356 590.6829
## Jan 2019       209.2025 133.5263 348.7313  94.31239 440.9304